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Despite the fact that the Schwarzschild and Kerr solutions for the Einstein equations, 
when written in standard Schwarzschild and Boyer-Lindquist coordinates, present coor- 
dinate singularities, all numerical studies of accretion flows onto collapsed objects have 
' been widely using them over the years. This approach introduces conceptual and prac- 

tical complications in places where a smooth solution should be guaranteed, i.e., at the 
gravitational radius. In the present paper, we propose an alternative way of solving 
the general relativistic hydrodynamic equations in background (fixed) black hole space- 
times. We identify classes of coordinates in which the (possibly rotating) black hole 
metric is free of coordinate singularities at the horizon, independent of time, and admits 
a spacelike decomposition. In the spherically symmetric, non-rotating case, we re-derive 
exact solutions for dust and perfect fluid accretion in Eddington-Finkelstein coordinates, 
and compare with numerical hydrodynamic integrations. We perform representative ax- 
isymmetric computations. These demonstrations suggest that the use of those coordinate 
systems carries significant improvements over the standard approach, especially for higher 
J - -^. ' dimensional studies. 
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Q\ . I. INTRODUCTION 

"o 

It is well known that the black hole solutions of the field equations of general relativity, with the 
| | Schwarzschild and Kerr solutions being astrophysically the most relevant, exhibit coordinate singularities 

when written in coordinates adapted to the exterior region. The very notion of a black hole was greatly 
clarified with the discovery, in the early sixties, of coordinate systems that remove those singularities, 
and indeed cover the whole spacetime |l| (for a general overview see 0], 0). The (rotating) black hole 
solution for the metric tensor, expressed in standard Boyer-Lindquist coordinates (t, r, 0, <j>) is singular at 
the roots of the equation A = r 2 — 2M r + a 2 = (where M is the mass and a the angular momentum 
aspect of the hole Q). In the spherically symmetric Schwarzschild case, the singularity at r — 2M is 
removable with the use of appropriate transformations of the radial and time coordinates. Slightly more 
complicated transformations, involving also the azimuthal coordinate, can remove the singularities at 
r = r± = M ± (M 2 — a 2 ) 1 / 2 of a rotating black hole. 

The location of the coordinate singularity coincides with the event horizon of the black hole. Asymptotic 
observers lose causal contact with events near the black hole at precisely this location. Hence, despite the 
singular appearance of the metric, in simulations of matter flows in the background gravitational field of 
a black hole, one is in principle allowed to consider the open interval extending from the event horizon to 
some "far zone" at a large, but finite, distance from the hole. This approach has been widely used over 
the years in the numerical simulations of flows around black holes || - jl2| . 

The blow-up of the metric components at the horizon has implications on the behavior of hydrodynam- 
ical quantities. The coordinate flow velocity becomes ultra-relativistic and reaches the speed of light at 
the horizon. As a consequence, the Lorentz factor becomes infinite causing any numerical code to crash. 
Placing the inner boundary close to the horizon, required for capturing the effects of the relativistic 
potential, introduces large gradients in all hydrodynamical variables. The steep radial behavior makes 
the task of numerical evolution with a reasonable degree of accuracy challenging. 

Besides practical considerations, the approach of evolving only the exterior domain lacks of a well 
defined location for the inner boundary: the computation cannot include the horizon surface, but must 
commence at a location "sufficiently close" to it. Physically, the influence of the horizon region on the 
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solution will progressively red-shift away the closer one gets to the horizon. The validity of a certain choice 
though, must be continuously reasserted, as new flows or effects are being investigated, with careful tests 
of convergence, as the inner boundary is progressively moved inwards. Such tests are complicated by the 
fact that, as mentioned in the previous paragraph, the solutions appear singular at the horizon and hence 
demand increasingly more resolution. 

Ameliorating those problems has motivated the use of a logarithmic radial coordinate (the so-called 
tortoise coordinate |l3fl). This technique relegates the event horizon to a infinite, negative, coordinate 
distance. An equidistant grid in the tortoise r* coordinate maps into an increasingly dense grid in the 
Schwarzschild r coordinate, with infinite density at the horizon. This approach has proven successful in the 
extensive semi-analytic studies of black hole perturbation theory, and recently also for the axisymmetric 
integration of curvature perturbations as an initial value problem jjflj. In wave systems, the ambiguity of 
the location of the inner boundary is addressed by the simple limiting form of the governing equations near 
the horizon. This is the well known fact that black holes act as finite potential barriers to electromagnetic 
and gravitational perturbations |l5| . 

The use of a tortoise coordinate does not bring similarly extensive benefits to the study of the hydro- 
dynamical equations. The issue of the inner boundary location is less transparent for those equations. 
The steepness of the solution and the "artificially" high coordinate velocities persist, although they are 
now more treatable due to the substantial increase in resolution. In three dimensional simulations using 
Cartesian coordinates, the tortoise technique is not possible at all. It has been shown recently that 
adaptive mesh refinement can indeed provide, even for 3D systems, the required resolution close to the 
event horizon. However, as alluded to above, these undesired pathologies can be eliminated in the case 
of fixed given black hole backgrounds with rather simple coordinate transformations. This reserves the 
power of adaptive mesh refinement for the more physically interesting features of the solution. 

There is considerable freedom in choosing coordinates regular at the horizon, which can be productively 
reduced by imposing criteria that can enhance their suitability for numerical applications. The obvious 
first criterion is of course the regularity of the metric, in particular at the horizon. A second condition 
is that the constant time surfaces are everywhere spacelike, as this is, currently, a pre-requisite for the 
implementation of modern numerical methods for relativistic hydrodynamic flows. An important third 
criterion is the time-independence of the metric components. This leads to constant in time coefficients in 
the equations and simplifies disentangling the true hydrodynamical evolution from coordinate effects in 
the black hole background. Interestingly, we will see that those conditions still do not fix the coordinate 
system uniquely. We show that the number of available choices greatly reduces, but is still infinite, in 
the rotating black hole case. We give several examples of such systems, which we collectively call horizon 
adapted coordinate systems. 

Such coordinate systems address the issues raised above in a straight-forward way: any radius inside 
the horizon is equally appropriate (in the idealized continuum limit) for the imposition of a boundary 
condition, as the domain is causally disconnected from the exterior. Importantly, the irrelevance of the 
inner boundary location will persist even after the inclusion of other possible local physical processes that 
may be considered in conjuction with the hydrodynamical flow, e.g., radiative processes. The coordinate 
velocities of the flow will be bounded at the horizon, as they represent projections of the (finite) fluid 
four velocity onto a regular coordinate system. Hence the hydrodynamical nature of the flow becomes 
considerably less demanding on the integration algorithm. Gradients in the solution for scalar variables 
such as pressure and enthalpy will of course persist. Those are physical and due to the curvature of the 
black hole, which requires a significant dynamic range for its resolution. 

The organization of this paper is as follows. In section II we introduce a class of horizon adapted 
coordinate systems for a non-rotating black hole. The rotating black hole and the more restricted class 
of coordinates available in that case are discussed in the Appendix. We outline the numerical hydrody- 
namical framework of our computations in section III. Exact solutions for spherical (Bondi) accretion are 
presented in section IV for both dust and perfect fluids. Section V describes the numerical results. In our 
numerical simulations we focus mainly on the spherically symmetric case, which captures the essential 
nature of the problem. Some axisymmetric computations are also briefly considered. The coordinate sys- 
tem on which we base our computations is the celebrated Eddington-Finkelstein coordinate system p7[ , 

Our main aim in this report is to show the functionality of this class of coordinate systems as frame- 
works for the integration of the equations of relativistic hydrodynamics in black hole spacetimes. In 
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three-dimensions (and rotating holes) computations of accretion onto black holes are likely to benefit 
significantly by the adoption of a horizon adapted coordinate system. 



II. HORIZON ADAPTED COORDINATE SYSTEMS 



A transformation of the Schwarzschild time t coordinate to a new "generic" coordinate t according to 

dt = dt±^(l-G(r)A(r)) 1/2 dr, (1) 

where A(r) = 1 — 2M/r, and G(r) an arbitrary function, brings the Schwarzschild metric into the form 

ds 2 = -A(r)dt 2 ± 2 (1 - G{r)A(r)) 1/2 dtdr + G{r)dr 2 + r 2 {d0 2 + sin 2 9d<t> 2 ). (2) 

The four-dimensional metric is next split following the {3 + 1} framework, as 

ds 2 = -{adtf + ll3 {dx l + I3 l di){dx° + ftdt) (3) 

and hence into the lapse function a = G(r) -1 / 2 , the shift vector /3, = (± (1 — G(r)A(r)) 1/ ' 2 , 0, 0) or 
/3 l = (3j, and the intrinsic spatial metric of the constant t hypersurfaces 7^ = diag(G(r), r 2 , r 2 sin 2 9). 

We proceed to establish the features we stated above as defining the horizon adapted coordinate 
systems. The spacelike nature of the foliation is evident for any positive definite function G(r), the 
normal to the hypersurface given byn^ = (G(r) 1 / 2 ,TG(r)- 1 / 2 (l-G(r)A(r)) 1 /2 j 0,0). The independence 
of all metric functions from the time coordinate t is also manifest. The regularity of the metric is also 
guaranteed if G(r) is bounded from above (we allow exceptions at the true curvature singularities). Hence 
a large class of horizon adapted coordinates is parametrized by the positive bounded functions G(r). The 
essential action of the coordinate transformation (|l]) is to introduce a non-diagonal g tr term (a radial 
shift vector). The shift covector will always have fi r = ±1 at the horizon and the time flow vector 
(#* = an^ + will be lightlike there. 

The choice of sign in the above formulae affects whether it is the past (-) or future (+) component of 
the event horizon that is regularized with the coordinate transformations. Some obvious choices for G(r) 
are G(r) = 1, G(r) = (1 + 2M/r), G{r) = (1 + 2M/r)(l + (2M/r) 2 ), etc. The second choice leads to the 
well known Eddington-Finkelstein (EF) coordinates which we will adopt for the numerical work in this 
paper. The coordinate transformation (]]]) and the relevant discussion of regularity is easily generalized 
for a rotating black hole. As shown in the Appendix, in that case additional conditions of regularity on 
an off-diagonal spatial metric component considerably reduce the class of horizon adapted coordinates. 
Further discussion of special cases is given there. Simulations of accretion flows onto rotating black holes 
in horizon adapted coordinates, more specifically in the Kerr-Schild coordinate system, will be presented 
elsewhere fl9| . We concentrate now on the EF form of the metric, regularizing the future component of 
the event horizon, which we will use for our hydrodynamical evolutions. With G(r) = (1 + 2M/r) the 
metric reads 

ds 2 = -(l- —) dP + —dtdr + (l + —) dr 2 + r 2 (d8 2 + sin 2 9d<j> 2 ), (4) 
\ T J r \ r J 

and correspondingly a = ( r+ ^ M ) 1 / 2 , f3 r = (2M/r) and 7^ = diag(l + 2M/r, r 2 , r 2 sin 2 9). 

The EF form of the metric represents an analytic extension of the Schwarzschild solution from the 
region 2M <r<ooto0<r<oo. It is among the algebraically simplest choices and generalizes 
to the case of a rotating black hole. In three dimensions it is naturally expressed in pseudo-Cartesian 
coordinates (Kerr-Shild form). An important coordinate property of the EF system is its adaptation to 
the null cone structure of the black hole: the ingoing radial null geodesies are described in a particularly 
simple way. 

Recent work clearly indicates that coordinate systems regular at the horizon, adapted to the local light 
cone structure, carry considerable advantages for the evolution of generic black holes. More specifically, 
the EF system and its generalization for dynamic spherically symmetric spacetimes has been used, for 
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the integration of a self gravitating scalar field in interaction with a black hole . In three dimensional 
integrations of the Einstein equations, the stable evolution of an isolated black hole has been an important 
target in the field of numerical relativity. The EF coordinates have been identified [glj as a point 
of departure for such integrations of the Cauchy problem. Recently, the Binary Black Hole Grand 
Challenge Alliance project demonstrated |22| the first advection of a black hole in a three-dimensional 
computational grid using the EF framework in boosted form. More remarkably, "eternal" evolutions have 
been achieved p3| using a formulation of the Einstein equations based on ingoing characteristic foliations. 

Those features of the EF family of coordinate systems do not appear particularly relevant for hydro- 
dynamical computations, since generic fluid flows are not tied to the spacetime null-cone structure in 
any essential way. It is apparent though, that in simulations accounting for the bidirectional coupling 
of matter to the spacetime geometry, computations of black hole interactions with matter would benefit 
significantly from this framework. This motivation breaks the degeneracy among coordinate systems for 
our purposes. 



III. RELATIVISTIC HYDRODYNAMICS AND NUMERICAL METHOD 

To demonstrate the feasibility of the procedure for numerical studies, we compute the hydrodynamic 
accretion of fluid flows onto black holes in the EF coordinate system. To do this, we solve the general 
relativistic hydrodynamic equations in these particular coordinates. Recently, Banyuls et al. pd| ] have 
written these equations as an explicit hyperbolic system of conservation laws (with sources) for a general 
metric written in {3 + 1} form. Here, we have to specialize those general expressions to the metric given 
by Eq. (Q). The evolved quantities are the relativistic densities of particle number, D, momentum (in 
the i direction), Si and energy, r. These quantities are defined in terms of a set of primitive variables 
(p, Vi,e) as D = pW, Si = phW 2 Vi and r = phW 2 — p — D. Here, p is the rest-mass density, p is the 
pressure, Vi is the 3- velocity of the fluid and h is the specific enthalpy, defined as h — 1 + e +p/p, with e 
being the specific internal energy. 

The contravariant components of the 3- velocity are defined as 

i u * P 

v l = — + -, (5) 
au l a 

where it' 1 is the fluid 4- velocity, a is the lapse function (a 2 — —l/g u ). Indices for spatial vectors (the shift 
vector and the 3- velocity) are raised and lowered using the spatial part of the metric, 7^. The quantity 
cm* is the Lorentz factor, W, which satisfies W = (1 — v 2 )~ x / 2 with v 2 — r )ijV t v^ . 

We have written a numerical code to solve the general relativistic hydrodynamic equations in the 
background spacetime of a non-rotating hole, using the EF coordinates. For this purpose we have taken 
advantage of the fact that the equations form a hyperbolic system and are explicitely written in conserva- 
tion form pd| ] . This allows for the use of advanced numerical methods based on approximate (linearized) 
Riemann solvers whose improved ability to handle relativistic (and ultrarelativistic) flows, high Mach 
number flows and discontinuous solutions (shock waves), have been recently investigated (see, e.g, pT| , 
JL2|| , p4| and references therein). The Riemann solver used (either Roe |29] or Marquina |2q] , [p7|) relies 
on the spectral decomposition of the Jacobian matrices of the multidimensional system of relativistic 
equations. The spatial accuracy of the code is improved by means of a monotonic piecewise-linear recon- 
struction |^s[| of proper rest-mass density, flow velocity and internal energy. Integration in time is done 
using a total-variation-diminishing Runge-Kutta scheme of high-order, developed by Shu and Osher p9| 
We refer the interested reader to Jll]] and |3(]] for further details on the system of equations and the 
numerical code. 



IV. EXACT SOLUTIONS 

The feasibility of the procedure is shown with a comparison of the numerical integration with exact 
solutions for spherically symmetric accretion of dust and perfect fluid. Here, we re-derive those exact 
solutions in the EF coordinate system. 
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Following Michel [[Uj we integrate the steady-state hydrodynamic equations to obtain 

V^gpu r = Ci, 



(0) 



~gpu r u t = C 2 , 



(7) 



where C\ and C2 are integration constants and g stands for the determinant of the four-metric, y/—g = 
a V7; with 7 = det (7^). Hence, Ut = C2/C1 = C3. The coordinate velocity v r is determined by C3 and 
u r through the condition g^u^u" = —1. 

For the dust case, the final expressions are simple algebraic functions. Fixing u t — C3 = —1 (the 
marginally bound case) the solution reads 



p(r) 



-Ci 



(8) 



v r {r) 



(9) 



which is valid in the domain < r < 00, i.e., also inside the horizon. The radial dependence of the rest 
of the hydrodynamical quantities can be directly computed from these two equations. 

Let us turn now to the perfect fluid case. Monotonic steady-state infall of perfect fluid onto a black 
hole must pass through the critical point of the so-called "wind" equation, obtained by differentiation of 
the integrated steady state equations (Eqs. (Q) and (Q)) with respect to p, u r and r, and elimination of 
the dp differential, 



du r 
u r 



V 



K) 2 



dr 
r 



2V 1 



M 



0. 



(10) 



where V 2 = ^On^ — ^ anc ^ u * = (" r ) 2 — 9 U ' Using the critical condition, that is, the assumption that 
both brackets in Eq. ( ^fj| ) vanish simultaneously, the relation (u r c ) 2 = M/2r c determines u r c and (wt) c , 
after the selection of the critical point radius r c . From u r we get Ut and hence we obtain V 2 at the 
critical radius, V 2 = (u r c /ut c ) 2 ■ On the other hand, considering a polytropic equation of state P — Kp r , 
with adiabatic index T and a critical density p c , we can infer the value of V 2 , which, in conjuction with 
V 2 = (u r c /u ti: ) 2 determines the polytropic index K. This complete determination of the fluid parameters 
at the critical point, combined with the steady-state conditions, fixes the integration constants Ci and C2, 
and hence uniquely identifies the extension of the solution away from the critical point. Contrary to the 
dust case, the solution does not have a closed form, as one must solve a non-linear algebraic equation for 
u r (r), from which one can get the rest of the variables. This is done using a Newton-Raphson iteration. 



V. RESULTS 



We now use the exact solutions derived in the previous section to check the accuracy of our numerical 
evolutions and to demonstrate the feasibility of our approach. For the numerical investigation we consider 
homogeneous initial data across the domain of integration. The initial data evolve in time, until a final 
steady-state solution is reached. We focus mainly on the spherically symmetric case, to which the exact 
solutions apply. At the end of the section we describe some axisymmetric computations. 

In spherical symmetry, the numerical domain extends from any given non-zero radius inside the horizon, 
fmim to some r max located in the physical universe (outside the horizon). In Figs. 1 and 2 we plot a 
comparison of the exact and numerical solutions for the dust and perfect fluid cases, respectively. We 
show the final steady-state for the spherical accretion problem. The solid lines represent the exact solution 
and the circles indicate the numerical one. In the particular computations considered here we have chosen 
Tmin = 0.5M and r max — 50M. We have used a non-uniform numerical grid of 200 zones (logarithmic 
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spacing). In both the dust and perfect fluid integrations we have only assumed the exact solution in 
the outermost zone. The integration has been performed until the final steady-state solution has been 
achieved in the whole domain. The perfect fluid equation of state is given by p = (T—l)pe, with T being 
the (constant) adiabatic exponent of the gas. We have chosen T = 4/3. 

The solutions are seen to be regular well inside the horizon at r = 2M. Nothing special happens at 
the horizon where all quantities are smooth. The steepness of the hydrodynamic quantities dominates 
the solution only near the real singularity. It is worth to stress the behavior of the velocity in this new 
coordinate system, which differs considerably from the Schwarzschild case. Now the total velocity does 
not approach the speed of light as it reaches the horizon, and although the maximum fluid speed is still 
at r = 2M, it is now considerably smaller, v = yj^ rr v r = 1/^3, in the dust case. For the perfect fluid 
case this value is even smaller, due to the pressure support of the gas. This is a most attractive numerical 
behavior, as it eliminates the divergence in the Lorentz factor at the horizon which would cause any 
relativistic fluid dynamics code to crash. As can be deduced from both figures we have found very good 
agreement with the exact solution. In Fig. (^|) we display the relative errors for density, pressure and 
total velocity in the case of perfect fluid accretion. 

The resolution benefits from the use of horizon adapted coordinates appear to be considerable, hence 
we attempt a rough quantitative assessment. We perform a sequence of comparison integrations in 
Schwarzschild and EF coordinates. We focus on the case of spherical dust accretion. We use a grid of 100 
zones and evolve the stationary solution up to a final time of 100M. The outer boundary of the domain is 
fixed at 50-M. The inner boundary of the EF domain is fixed at 1.5M. For the Schwarzschild computation 
we use a tortoise radial coordinate which provides very high resolution near the horizon. We have found 
that the solution in Schwarzschild coordinates is maintained stationary, provided r min > 2.1M. For 
Train = 2.05M and a maximum resolution of 10~ 3 M, the accretion departs from the exact initial solution 
at the zones closer to the horizon. For r m i n = 2.01 and same resolution, the code crashes. In this case, a 
grid of 200 points does not preserve the initial solution either (even with a resolution of 10 _4 M). At least 
400 zones are required (with a maximum Ar 10 -5 ) to achieve the correct accretion pattern. On the 
other hand, using EF coordinates, we can resolve the problem with an (unequally spaced) grid of 50 zones. 
This gives, broadly speaking, a factor of 8 benefit in resolution when using horizon adapted coordinate 
systems. Computations that may require small values of r m i n will further increase this number, the same 
being true for perfect fluid computations. In this latter case the computation in standard coordinates is 
more challenging, as the pressure (or the internal energy) can easily become unstable near the horizon 
and even become negative. Finally, three dimensional Cartesian computations do not allow the use of 
tortoise coordinate based grids. Hence a conservative estimate introduces an 8 3 gain factor in the size of 
tractable problems in three dimensional integrations. 

We now briefly describe axisymmetric computations. We study the so-called Bondi-Hoyle accretion, 
that is, hydrodynamic accretion onto a moving black hole, in the EF coordinate system (for a general 
overview of such flows see, e.g., [ fL2f and references therein). The asymptotic values of fluid velocity, Voc, 
sound speed, c Soo and adiabatic exponent, T, together with the assumption of an initially homogeneous 
distribution suffice to define the initial values of all hydrodynamical quantities. As in the spherical case, 
these data evolve towards a final steady-state solution. We use a grid of 200 x 100 zones in r and 0, 
respectively, with r m i n = 1.5M, r max — 100AZ and < 8 < it. The radial grid is logarithmically spaced. 
We evolve a model with c Soo =0.1, and = 0.5. Hence, it is a supersonic model with Mach number 
5. We plot in Fig. 3 the stationary solution for T = 4/3 (top), 5/3 (middle) and 2 (bottom). We plot 
iso-contours of the logarithm of the density. The most prominent feature of the solution is the stationary 
conical shock. The inner circle indicates the position of the horizon. Again, as in the spherical case, 
the solution smoothly straddles the horizon at 2M. Compared to the evolutions presented in Jl^] in 
Schwarzschild coordinates, we now use less resolution to get to the same final solution. Highly supersonic 
models with large values of T, in particular T = 5/3, were very hard to evolve with Schwarzschild 
coordinates and the T = 2 case was impossible to evolve. Its computation is now straightforward. 

VI. CONCLUSIONS 

We have presented a family of horizon adapted coordinate systems for the numerical study of accretion 
flows around black holes. In the rotating case we identify a discrete but infinite family, of which the 
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first simple members are well known stationary coordinate systems. In the non-rotating case the freedom 
in building the regular stationary foliation is quantified by (at least) the space of bounded positive 
functions of one variable. We have shown how these systems allow for a better numerical treatment of 
accretion scenarios. Existing numerical studies of accretion flows onto black holes have been performed 
in the original, singular systems, i.e., Schwarzschild coordinates for the Schwarzschild solution and Boyer- 
Lindquist coordinates for the rotating (Kerr) solution. Although it is possible to solve the problem in 
these pathological coordinates, one is introducing artificial complexity, being forced to use very high 
resolution to deal with the unphysically large gradients that develop in the vicinity of the horizon. This 
may prevent the accurate solution of three dimensional problems. At the same time, the ambiguity 
regarding the position of the inner boundary of the domain (which should be the horizon) introduces a 
convergence criterion that must be enforced at all times if the solutions are to be trusted. 

We have focussed on the particular case of the Eddington-Finkelstein form of the Schwarzschild metric. 
The general relativistic hydrodynamic equations are now regular at the horizon, which permits an accurate 
description of accretion flows. We have demonstrated the feasibility of this approach with the numerical 
study of the spherical accretion (Bondi accretion) of dust and perfect fluid, and the comparison with the 
exact solutions which we re-derived in this coordinate system. We have also shown the functionality of 
the new coordinates in axisymmctric computations of relativistic Bondi-Hoylc accretion flows. 

In a forthcoming paper we plan to extend this approach to the rotating case, considering horizon 
adapted coordinate systems to study accretion flows in stationary Kerr spacetimes. Three dimensional 
accretion flows onto black holes are interesting both from an astrophysical and geometrical point of view, 
as they are thought to correspond to observable electromagnetic emission, and hence may help map 
the relativistic rotating black hole potential. The framework proposed here will help detailed numerical 
studies of such systems in the near future. 
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APPENDIX: 



In this appendix we present stationary coordinate systems for the Kerr spacetime that are regular at 
and outside the horizon r+, i.e., the larger of the roots of the equation A = r 2 — 2Mr + a 2 = 0. In the 
Boyer-Lindquist form for the metric tensor the j rr component is singular at r + , 



ds 2 



A 



— [dt — a 



sin 



2 sin t) r , 2 2s , , „-, 2 
H — [(r 2 + a 2 )d(j) - adt\ 



P 2 



+ ^dr 2 +p 2 d0 2 , 



(Al) 



with p 2 



■ a 2 cos 2 9. We now introduce a coordinate transformation defined by 



d(f> = d(f>— -^ e ( r 7 0)dr, 



di = dt- e(r, 6) 



r 2 + a 2 



p 2 - 2Mr 



2Mr 



G(r,6) 



dr, 



(A2) 
(A3) 



where e and G are, at present, arbitrary functions. This form of the transformation is chosen in order 
to isolate the singular behavior of the spatial metric. With the above ansatz, the metric becomes, in the 
new coordinates (t, r, 8, <j>), 
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ds 2 = - 1 - 



2Mr\ 



dt 



72 



AMar 



sin 2 Odtdcj) 



(A4) 



-2c ,/l -II- | GdWr 



1 -e z 



p 2 + G 



+2e 



a sm 



2 fl„2 



p 2 - 2Mr 



+p 2 d6 2 + sin 2 6> 



a z sin z 6» 1 



2Mr\ 
—J 



dr 2 



dcfidr 



By inspection of the spatial metric components we see that the condition of regularity at the A = 
surface demands that e 2 — 1. This fixes the function e up to a sign. We proceed to further examine the 
regularity of the metric element. All components of the metric, except j r ^, are regular, except at the 
points p 2 — 0, which belong to the true singularity. The component 7^ is singular at the ergosphere 

(p 2 — 2Mr = 0) which is outside the r + horizon, and coincides with it at the poles 9 — 0, ir. 

Before we ameliorate this pathology we note that in the case of no rotation (a = 0), the off-diagonal 
term drops, hence as we have already seen in section there is no restriction on the choice of the 
function G. This special behavior of the Schwarzschild black hole allows selecting G with the target of 
simplifying the spatial 3-metric. An interesting coordinate system, is obtained with the simple choice 
G(r) = 1, which brings the metric into the form 



/ 2M\ ~ 2M ~ , , 
ds 2 = - ( 1 jdt 2 + 2d dtdr + dr 2 



and correspondingly a = 1, /3 r = (2M/r) 1 ^ 2 and 7^ = diag(l, r 2 , r 2 sin 2 6). This is the Lemaitre 
coordinate form of the Schwarzschild solution |32|]. The trivial form of the lapse and the three-metric 
imply that all of the curvature information of the black-hole spacetime is encoded in the way the r, 9, (f> 
coordinates must propagate from time slice to time slice (choice of shift vector), in order to preserve 
a flat 3-space, and a constant rate of clock ticking. It is quite remarkable that such a shift vector is 
a simple, time- independent, function of radius, which, incidentally, is equal to the Newtonian free fall 
velocity as measured locally by static observers. This form of the metric element was the starting point 
for a program of constructing black hole initial data p3| . Unfortunately, this simple state of affairs does 
not appear generalizable to the rotating case. 

Returning to the general case, by further exercising the freedom in choosing the function G(r, 9) we 
can ensure that the off-diagonal metric component is regular. 
Substituting for G via 



r 2 {d9 2 



(A5) 



/ 2Mr\ , 



(2Mr\ 



we see that the class of regular systems is parametrized by the functions 

fe 

{^) ' 

where k is a non-negative integer. The general admissible form for the function G(r, 9) is 

2fc-l 



G fe = 1 + 



2Mr 



/2MrV 



(A6) 



(A7) 



(A8) 



for non-zero fc, and 



Go = 0. 



(A9) 
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With those substitutions the metric becomes 



ds 2 = - 1 - 



2Mr\ , ,, 

— r 



4Mar 



sin 



a 2 sin 2 



sin 2 9did<j) 

2Mr 



(A10) 



1 



-2eZ k dtdr + G k dr 2 + 2easin 2 9W k dcj)dr, 



where 



2Mr 

T 2- 



(2Mr 



(All) 



and all terms are regular except at the true curvature singularity. 

The case k = corresponds to the ingoing and outgoing Eddington-Finkelstein coordinates (e = — 1 
and e = 1 respectively). In this case the coefficient of the 7 rr vanishes. Despite a common misconception, 
the foliation is still spacelike H] for non-zero a. The null character of the foliation in the limit a = 0, 
though, implies that this coordinate system is not appropriate for 3 + 1 based studies. 

The case k = 1 generates the Kerr-Schild coordinate system, expressed in spherical coordinates via the 
relations 101 



x + iy — (r + ia) sin 9e 1 ^ , 
z = r cos 9. 



(A12) 
(A13) 



It is the algebraically simplest member of the family and reduces to the EF system in the a — case. 

The case k = 2 introduces a time coordinate that satisfies a harmonic condition ]35| ], i.e., Dt = 0. Some 
recent hyperbolic formulations of the Einstein equations are based on harmonic slicings, hence the black 
hole metric expressed in such coordinate systems can be used effectively as a numerical code test-bed Q . 

The higher order members of the family, i.e., k = 3,4, .. . may also have interesting geometrical or 
differential properties. 
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FIG. 1. Exact (solid line) versus numerical (circles) solution for spherical accretion of dust, after reaching 
the steady-state. The top curve shows the total velocity and the bottom one shows the rest-mass density, as a 
function of the radial coordinate. The total domain extends to 50M. We focus on the first 15M. The value of 
Ci in Eq. (^) was chosen to be —0.195. The numerical solution agrees well with the exact one. The dotted line 
marks the position of the horizon. 
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FIG. 2. Exact (solid line) versus numerical (circles) solution for spherical accretion of a perfect fluid, after 
reaching the steady-state. The top curve shows the pressure, the middle one shows the total velocity and the 
bottom one shows the rest-mass density, as a function of the radial coordinate. As in Fig. 1, the total domain 
extends to 50M and we only plot the first 15M. The numerical and exact solutions also match very well. The 
dotted line marks the position of the horizon. The critical point values used to construct this particular exact 
solution are r c — 400M and p c = 1CF 2 . 
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FIG. 3. Relative errors for density, pressure and total velocity in the case of perfect fluid accretion, as a function 
of the radial coordinate r. The result is for 200 grid points between 0.5 and 50M. Only the first 15M are shown. 
The overall good accuracy is apparent, the especially small error in total velocity being most striking. 
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FIG. 4. Grey-scale iso-contour plots of the logarithm of the rest-mass density for the axisymmetric relativistic 
Bondi-Hoyle accretion problem. The final stationary solution is plotted. From top to bottom we display the 
solution for T = 4/3, 5/3 and 2, respectively. The left column shows a close-up view of the central region, which is 
enlarged in the right column. The inner circle shows the location of the horizon. The shock cone is well resolved 
even at large distances from the hole. Its innermost location differs in the three cases. For larger values of F it 
moves to the front part of the hole, due to the higher pressure values at its rear part. 
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